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Abstract 

Nonlinear dynamical systems are sometimes under the influence of random 
fluctuations. It is desirable to examine possible bifurcations for stochastic 
dynamical systems when a parameter varies. 

A computational analysis is conducted to investigate bifurcations of a 
simple dynamical system under non-Gaussian a— stable Levy motions, by 
examining the changes in stationary probability density functions for the 
solution orbits of this stochastic system. The stationary probability density 
functions are obtained by numerically solving a nonlocal Fokker-Planck equa- 
tion. This allows numerically investigating phenomenological bifurcation, or 
P-bifurcation, for stochastic differential equations with non-Gaussian Levy 
noises. 
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1. Motivation 

The dynamical behaviors for a dynamical system depending on a pa- 
rameter may change when this parameter varies. This so called bifurcation 
phenomenon has been observed in many deterministic systems It also 
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occurs in stochastic systems with Gaussian noises |3j. It is desirable to exam- 
ine possible bifurcation phenomena for stochastic systems with non-Gaussian 
noises. 

Levy motions L t are a class of stochastic processes that have independent 
and stationary increments. They are usually non-Gaussian processes. The 
well-known Brownian motion B t is a special case which has additional prop- 
erties: (i) Almost every sample path of the Brownian motion is continuous 
in time in the usual sense and (ii) Brownian motion's increments have Gaus- 
sian distribution. Random fluctuations in complex systems in engineering 
and science are often non-Gaussian. For instance, it has been argued that 



diffusion by geophysical turbulence (17| corresponds, loosely speaking, to a 



series of "pauses" , when the particle is trapped by a coherent structure, and 
"flights" or "jumps" or other extreme events, when the particle moves in the 
jet flow. Paleoclimatic data Q also indicates such irregular processes. 

SDEs perturbed by non-Gaussian Levy noises have attracted much at- 
tention recently 0, Oil]. SDEs perturbed by Levy motion generate stochastic 
flows 10, [2| , or random dynamical systems (cocycles), under certain condi- 



tions. 

Let us consider a deterministic differential equation perturbed by a non- 
Gaussian Levy motion, i.e., consider a stochastic differential equation (SDE) 

dX t = f{b,X t )dt + edL?, (1) 

where b e R, e > and a G (0, 2) are real parameters; and Lf is a a— stable 
symmetric Levy motion defined in a probability space (Q, J 7 , P). In this 
paper, we consider a numerical approach for understanding how the dynamic 
behaviors change when parameters vary, for a special case / = bX t — Xf. Note 
that x = bx — x 3 is a primary dynamical model exhibiting the deterministic 
pitchfork bifurcation [8[. 

In section 2, we briefly review some basic concepts for Levy motions. In 
section 3, we first discuss the nonlocal Fokker-Planck equations for SDEs with 
Levy motions, then present a numerical approach in computing stationary 
probability densities for the solution processes for (P) above, and further 
discuss how the stationary probability densities change when the parameters 
vary (i.e., phenomenological bifurcation or P-bifurcation). 



2. Levy motions and generators 

Let us briefly review basic facts about Levy motions. A scalar Levy 
motion is characterized by a drift parameter 9, a variance (or diffusion) pa- 
rameter d > and a non- negative Borel measure u, defined on (R, S(R)) and 
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concentrated on R \ {0}, which satisfies 

/ (y 2 A 1) u(dy) < oo, (2) 



or equivalently 

——2 u(dy) < oo. (3) 
1 + y 2 

This measure v is the so called Levy jump measure of Levy motion L(t). We 
also call (9, d, v) the generating triplet. 

Let L t be a Levy process with the generating triplet (9, d, v). It is known 
that a scalar Levy motion is completely determined by the Levy-Khintchine 
formula (See [I EH]). This says that for any one-dimensional Levy process 



L t , there exists a 9 e R, d > and a measure z/ such that its characteristic 
function is 

Ee iXLt = exp{i9\t - dt— + t / (e lXy - 1 - i\yI {M<1} )v(dy)}, (4) 

2 JK.\m 



where Is is the indicator function of the set S, defined as follows: 

is(y) = 



if y e S; 
My iS. 



The generator A of the process L t is the same as the infinitesimal gener- 
ator since Levy process has independent and stationary increments. Hence 
A is defined as Aip = lim^o where Pt<p(x) = E x tp(L t ) and ip is any 

function belonging; to the domain of the operator A. Recall the generator A 
for Lt is (See QQ) 



Acp{x) = 9<p'(x)+-(p"{x)+ / [<p(x+y)-(p{x)-I M<1} y<p'(x)] v{dy). (5) 

z JR\{0} 

In this paper, we consider a special Levy process, the symmetric a— stable 
Levy motion Lf, with drift 9 = 0, diffusion d = and the jump measure 
v a (dy) = jjirz- The corresponding generator is A a ip(x) = f R ^ {0} [ip(x + y) - 

- h\y\<±} yf'( x )]^a(dy). 

In the next section, we consider bifurcation of the equation ([T|) when b 
and a vary, by numerically investigating stationary probability density func- 
tion for the solution of ([T]). We take the drift coefficient f(b,x) = bx — x 3 , 
corresponding to the well-known pitchfork bifurcation in the deterministic 
case (when e = 0). Note that simulations for solution paths were conduced 
for SDEs with a— stable Levy noises in while we examine bifurcation phe- 
nomena by computing stationary probability density functions for solutions 
in the present paper. 
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3. Bifurcation under additive Levy noises 

We now consider possible bifurcations for the SDE 

dX t = f(b,X t )dt + edL?. (6) 
The generator A for the solution process X t in ([6]) is 

Aip = f(b,x)ip'(x)+e [ip(x + y)-ip(x)-I { \y l<1} y(p'(x)]iy a {dy). (7) 



The adjoint operator for A in the Hibert space L 2 (R), with the usual scalar 
product, is then [l[ 



A*<p = -[/(&> x )f( x )\ + e / [tp(x + y)- tp(x) - I{\ y \<i} y<p'(x)] v a {dy). (8 

Jr\{o} 

Consequently, the Fokker-Planck equation for the probability density func- 
tion p(x, t) for the solution process X t in (JTJ) is 

d t p = -d x [f{b,x)p{x,t)} 

+ el [p(x + y,t) -p(x,t) - I{\ y \<i} y d x p(x,t)]\ vjdy). (9) 



The stationary solutions p(x) of the above Fokker-Planck equation define 
some invariant measures for the equation ([T]): for real parameter b and Levy 
parameter a G (0,2). Namely, the stationary probability density function 
p(x) satisfies 

- [f{b, x)p(x)]' + e / [p(x + y)- p(x) - I { \ y \ <1} yp'(x)] —— = 0, (10) 



p(x) > 0, / p(x)dx = 1. (11) 
Jr 

3.1. Deterministic pitchfork bifurcation 

The differential equation it = bxt — x^ has a stable fixed point at x = 
for b < 0, and two additional stable fixed points at x — ±y/b for b > (See 
j8[). It undergoes a so-called pitchfork bifurcation at b = 0. 



3.2. Bifurcation under Brownian motion 

We first recall a bifurcation under Brownian motion (See [sj or 0, Page 
475]), i.e., in the case when the Levy motion in the above equation ()6]) is 
replaced by a Brownian motion W t . For dX t = {bX t — X^)dt + adW t , b G 
R, a 7^ 0, there exists only a unique stationary measure with density p^ = 
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N^p exp[^2 (bx 2 — where is a normalization constant. For any give 
noise intensity a 7^ 0, the density is unimodal for b < 0, but bimodal for 
b > (and the plateau for p(x) occurs at x\ = y/b and X2 = — Vb). Hence 
the family ( P(b,a) beR ) undergoes a bifurcation at b = for each a ^ 0. 
This is a kind of phenomenological bifurcation or P-bifurcation Q, in which 
a Brownian motion leads to a different bifurcation than its deterministic 
counterpart. 

3.3. Bifurcation under a— stable Levy motion 

Now we consider the bifurcation under a— stable Levy motion. For dX t = 
(bX t — X^)dt + edL",b G R, a G (0,2), a stationary measure has density 
p(x) = Pb,e,a( x ) satisfying the following steady Fokker-Planck equation 

- [(bx - x 3 )p(x)]' + e / [p(x + y)- p(x) - I {lyl<1} yp'(x)] --— = 0, (12) 

</r\{o} \y\ 



p(x) > 0, / p(x)dx = 1.(13) 
Jr 

Unlike the Brownian case in §3.21 above, we do not have the exact analyt- 
ical solution for the equation (fl2l) . In order to detect possible bifurcations, 
we instead numerically simulate this integro-differential equation on the in- 
terval (—1,1), with / > large enough and with the homogeneous Dirichlet 
boundary conditions. This integro-differential equation contains both an dif- 
ferential part — [f(b, x)p(x)} and an integral part J R \^[p(x + y) — p(x) + 

^■{\y\<^-}yP '( x )] 1 • ^ e use a fi m t e difference scheme on the differential 



part and the trapezoid rule in the integral part [l3[; see also [5]. 

We conduct numerical simulations for various b 6 R, a 6 (0, 2) and e G 
(0, 1). Although b and e may be any real number in our numerical approach, 
here we limit them b to be in a bounded interval in this paper. All figures 
are in color in the online version of this paper. Different colors are used to 
distinguish cases with various parameter values. In the following only some 
selected figures are shown to illustrate our results. 

3.3.1. Varying the parameters b and a 

Figures [1] and |2] show the stationary probability density function p(x) 
for e = 0.1 and e = 0.9, respectively. Here we only show several cases for 
b = -5, -1, 0, 1 and a = 0.1, 0.4, 0.7, 1.0, 1.3, 1.6, 1.9, 1.999, as examples. 

The probability density function p(x) evolves from bimodal to unimodal, 
and then further changes to the flatter kurtosis shape for every fixed pa- 
rameter b in interval(— 10, 0), as a value increases. The bifurcation occurs 
only for b < 0, at some bifurcation value a* = a(b, e), in our computational 
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Figure 1: Stationary probability density p(x) for e = 0.1 



range. Note also that the p(x) is bimodal for b > 0. This phenomenon is 
more evident when e is larger in the range (0, 1); see Figure 2. 

For fixed b value, the probability density function p(x) evolves from lower 
kurtosis to the higher one, and then changes to lower again, as a value 
increases. 

3.3.2. Varying the parameter b 

When a is approximately within in the interval (0.4, 1.6), the stationary 
density p(x) becomes very spiky, and this is more evident when the magnitude 
of b is large; see Figure |3j 

3. 3. 3. Varying the parameter a 

When a is approximately within in the interval (1.9,2), the stationary 
density p(x) becomes very flatter, and this is more evident when the value 
of a is close to 2; see Figure HJsubfigures (c) and (d) in both figures). 

3.3.4- Impact of noise intensity e 

When the positive noise intensity e is increased, we observe that the 
stationary density p(x) becomes flatter (or less spiky) for fixed parameters b 
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Figure 2: Stationary probability density p(x) for e = 0.9 
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Figure 3: Stationary probability density function p(x) for e = 0.1. (a) Corresponding to 
a = 1 and b = —30; (b) Corresponding to a = 1 and b = —50; (c) Corresponding to 
a = 1.5 and b = 20; (d) Corresponding to a ~ 1.5 and b — 50 




Figure 4: Stationary probability density function p(x) for e = 0.1 and a = 1.999. (a) 
Corresponding to b = —30; (b) Corresponding to b = —50; (c) Corresponding to b = 20; 
(d) Corresponding to b = 30 
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Figure 5: Stationary probability density p(x): a = 0.3, b = —10 for e = 0.1, e = 0.9 and 

e = 3 
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Figure 6: Stationary probability density p(x): a = 0.7, b = —10 for e = 0.1, e = 0.9 and 

e = 3 
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Figure 7: Stationary probability density p(x): a = 1.6, b = 5 for e = 0.1, e = 0.9 and 
e = 3 
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Figure 8: Stationary probability density p(x): a = 1.5, b = —8 for e = 0.05 and e = 0.1 
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and a; see FiguresOlEJ Thus the nonlocal integral term, due to Levy jumps in 
the random forcing, in the steady Fokker-Planck equation (1121) has a certain 
damping or diffusive effect. 

3.3.5. A remark: Brownian motion vs. a— stable Levy motion 

When a = 2, the corresponding Levy motion is the Brownian motion. 
We take a close to 2, the stationary density p(x) are computed for various b 
values; see Figures [3] and H] (subfigures (c) and (d) in both figures), and also 
Figure [3 Here we also observe that p(x) is bimodal for b > 0, similar to the 
Brownian motion case in §3.21 
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